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In thxa report our main concern u linear 
prediction as veil as spectral estination in two 
dimensions. She * 78 X 1008 models applicable for tvo- 
dimensional prediction ere £p.ven. Both coinreational 
and ISHSUH modd for 'two dimensional spectral estima- 
tion arc discussed. Xhe result includes many numenoal 
examples 9 a de-bailed bibliography and a listing of 
related oonputer routines. 



OHAPSER 1 


UnERODOOTIOH 


With the advent of digiteJ. oomputere and digital lOs, 

digital techniques for solving engineering problems have 

At 

become increasing popular. Using the digital approaoh, the 
designer no longer has to be oonoemod with realisabili'lT' 
constraints of analog devices. Digital eifstems arswconoerned 
with discrete time inputs and their processing to achieve a 
desired obsective. Before the highspeed digital computers 
were developed discrete time signal processing was limited 
to relatively simple algorithms. Advancements in digital 
comiuter technology revealed enormous potential of computers 
and digital signal processing has become a very Important 
part of many areas of science and technology. 

One of the most important areas of applioatlon, for 
digital signal processing techniques is the estimation of 
power spectrum of a time or space senes. OJhe need for 
power spectrum arises In a variety of contexts including the 
measurement of noise spectra for the design of optimal linear 
filters, the detection of narrow band signals in wide bend 
noise, the estimation of parameters of a linear system by 
using a noise excitation etc. 

Ihe present work consists of two parts, namely: 

(1) two dimensional prediction, (2) two dimensional spectral 
analysis. 
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She Bpeotral analysis of multi-di'oensional sigtials 
plays a very important role In areas, suota as, image prooessing, 
radio astronomical appliontions, geophysics and solid state 
pbysios* In radio astronomy itself it has many applloatlons: 
apertnre synthesis, interplan etaxy scintillation analysis and 
radar haok scatter studies. 

Sphere are basically two conventional methods of 
power spectrum analysis, namely the autocorrelation or the 
Blaokman-Iukeiy method and the "periodogram" method. Both these 
methods, which are nonparometrio , make certain assumptions 
about the data outside the observation interval and their 
resolution performance is directly dependent on the length 
of the observation Interval. Some of these limitations 
prompted the search for alternative schemes and one of the 
most promising one to cone by m recent years is the so-called 
"maximum entropy method". She method is oonsistant with the 
statistical information av*«ilable about the data and is 
"maximally non-committal" about tho data outside the observa- 
tion interval. The method has shown drama tic results in Iftie 
resolution capability on short data lengths, furthermore, 
the method yields results which are shown to be equivalent 
of fitting an optimum autoregressive or prediction model 
to the data under consideration. 

The present report, therefore, begins with a 
discussion of autoregressive or prediction models for two 
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dimeasioncG. data . Chapter 2 deals with various types of 
predictive filtering. In each case the "normal equations" 
are derived using Minimum Mean Square Error Criterion 
(MMSE). Justice [18] has given an algorithm, (based on the 
theory of Sseigo polynomials), whloh Is analogous to the 
Levinson's algorithm in one dimension. !Ehis has been 
disoussed in somo detail. Chapter 3 deals with conventional 
methods of spectral estimation. Chapter 4 deals with the 
maximum entropy method in 2-dimensions. Ihis chapter shows 
how Wiener filter theory can bo extended to two-dimensional 
oase and used for power spectral estimations. !Che concluding 
chapter presents and disousses some of the computational 
results. She computer programmes are listed in Appendix. 



OmSTES. 2 

LIinB&R PBEHOmON IN WO DIMENSIOMS 


Xhp matheoatioal analysis of gensral dynamo ey stems 
bas been an area of masor oonoem slnoe the beginning of this 
oentuxy. With the advent of digital oomputerSf ibis problem 
has been pursued with greatest vigour. She aanalysls of 
outputs of dynamo ay stems has been of fundamental importanoe 
in the fields of eommuni cations, statistics and eoonomos. 

She ' time series* analysis as it Is called is used in maoy 
areas of oominuni cations sudh as: speech and picture trans* 
mission, automatic speech recognition and data compression. 

An important oonoem of time series analysis is the 
estimation of power spectra, cross spectra, and coherence 
functions. An Important technique In many of these studies 
IS * Linear Prediction*. In linear prediction, we try to fit 
a linear model for the process under study, based on the 
statxiBtios of the observed teta. One suoh model Is the Auto 
Regressive Moving Average (ARKA) model. In ARKA model, the 
signal Sq is considered to be the output of a eystem with some 
unknown input u^, such that the following relation holds good. 

®n ■ “ jgJi ®k ®n-k * ^ ^n-1' ^o “ (2.0.1) 

where a|g, lilk;! p, b-j^, l£. q, and the gain 0 are the 

parameters of the hypothised system. Equation (2.0.1) says 
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that the output at waj instant of time is a linear 
oombination of past outputs and present and past inputs* 

Shat IS the signal is prediotable from 'ttie linear 
combination of past outputs and inputs. Henoe the name 
* Linear prediction^. Having fitted an ABHUL model for the 
process, we oan use the parameters of this AEMA model for 
different applications such as forecasting, data compression 
and power spectral estimation. 

She philosopb7 of linear prediction oan be extended 
to two or more dimensional data. Here instead of a signal 
with one independent variable n, we have a signal Sj^^^ 
with two spatial variables k and 1. ^ making use of 

suitable strategies, it is possible to develop a set of 
equations which arc similar to that in one dmoisional case. 
These equations oan be solved recursively to obtain the 
parameters of the model. 

2. 1 PHEDIOTIQH HOHSL PQR SOAlTiED lUTA; 

In this and the following sections, we discuss 
various models suitable for two dimensional prediction* The 
one under oonsideration is suitable for scanned two dimensional 
data such as pictures [11]. Let us imagine that the given 
picture IS scanned into H sequences which are stationary In 
nature (Pig.1). Let Z(i,n) denote the nih sample of the ith 
sequence* Let the cross correlation between two sequenoes 
be defined as 

1(1, l+Wl 


( 2 . 1 , 1 ) 



m, 

-H 

►T 


M 

•i 



Fig.l Sehonatlo Arxangeffleat for Prediction Stxcs.teg7 
uaed in Section 2>1» 
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!Che prediction strategy is 
R N 

l(i,n) ml Jo (i,3fk) I(3fn-fc) (2.1.2) 

to-1 3-1 


where a (iy 3 ,i£) are constants to be suitably chosen. ITow 
we choose the coefficients such that lAie mean square error 
between the actual samxd.e X(i,n) and the predicted values 

A 

Z(i,n) IS minimum. 


Xhus we ire required to minimise the mean square 
error given by 


e*- E [(iCi.n) - i(i,n)} 

R R 

B l<X(i,n) -II o(i,3,k) 1(3, n-k)) ^ 1 

Itel 3=1 

(2.1.3) 


Setting 


3 * T 


3* r 


3a (1,3 


B 0, we have 


R R 


- 2E C-I( 3 , n-k) {X(i,n) - J X o(l,itt,q)I(m,n-q)}l 

3a (l, 3 * 1 c) q-1 mtl 

- 2 a(i,m,q)5jnj(q-k) 1 2 0 

Xhus a(i,j,k)*s are the solutions of the following equation. 

R.,(k) -if a(i,m,q) 

31 q-1 m-1 


(2.1.4) 


iEhus there are unknowns with equations and therefore 


these equatiohs Odn be solved. Uhls formulation represents 
the two dimensional generalisation of the prediction problem 



7 


Let us try to pat this in the matrix form. Xoward this 
end, let 

Hjj(k) - S(i,u) 

o B (!,■») 

- T(t,u) 

iMhcre u is any unique one to one mapping of the ordered pairs 
(Stk) and ▼ is the same mapping of the ordered pairs (m,q). 
Vor exeunple, let 

XL - (3-1 )H + fc 

▼ « (m-1)R + q 

I 

She prediction equation oan then he written as 


[S(i,u)] « [!(▼, tt)I [0(1, ▼) 1 (2.1.5) 

2 

where X(u,v) is a matrix with IT suhmatrioes as shown below. 
[IL|^(k-q)] [R2j(k-q)] ... [R|f^(k-q)] 

[R22(^s-4)] ••• 


T(u,v)- 


U [R^jj(k-q)l [R2jf(k>q)] ... tR|pj(k-q) ij 


Xhus I(u,v) consists of submatnoes which are loepHts 
matrioes. Xbus the equation (2.1.5) oan be solved by using 
algorithms for inverting loeplits matrioes [02] , 
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2.2 VIMSR FinESRlBG MOSEL: 

In this section, our main concern is to develop a two 
dimensional Wiener filter. She Wmer filter should be such 
that given an input arraj b ■ must be able to obtain 

the desired output array o m {Oj^} as the convolution of the 
input array b and the filter array a « ^^le formulation 

of the filter equations is as follows C18] . 


We are given arrays b 

satisfying 
08 00 

I I ^ 

k-o 1-0 


{b^^} and o - ( Oj ^} , k >0, 1 iP 

00 00 « 

I I lo^l • 

k— o 1—0 


Our problem is to find a (M+l ) b7 (M+1 ) array »a* for whibh 
it 18 true that the convolution b « a is the best least squsires 
approximation to the array Making use of the fact that Ihe 
spaoe of square summable sequences is a Hilbert space, we oan 
derive a set of equation whose solutions yields the unknown 
array a. 

Lefino arrays u^_ as the translates of the b array 
i.e. U|j2l(1i3) - £i-1c,D- 1 negative subsonpts)- 

We define an inner product by 

■w* 

< a,b siyi b^^ 

and norm in the usual way» 

|a||^-<a,a> -EE I Uiri I 
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Oar pro'blem is Ibeii to minimi se the norm 


n 

II i 

i 


H 

J. 

O l£mO 


®kl ® I 


Ibat is, wo must find the best approximation to o from the 
suhspaoe spanned hy the linearly independent amreys • 

It IS well known that the solution satisfies the equation 
N M 

"J J aicl"kl‘®' for OiB iH, 

KhO laO 

0:lt <H 

That IS, the difference array is orthogonal to the indioated 
suhspaoe. Ve may rewrite the equations in the foxn 


ST N 

®kl ♦B-k.W. “ «Bt OS.^, 0 SUM 

( 2 . 2 . 1 ) 

!Chese are "Normal equations". 

The array g^^ «< o, Ug.^> Is the "oross-oorrelation” 
of the desired array o with the input array !Che array 

the "auto-corr elation" of the input array. 
It oan he easily verified that a unique solution { a^} exists 
for any aooeptahle choice of the array o. Also the definition 
of the oooffioient implies that 

♦kl " ♦-k,-l 

Let us denote the row vectors of a ■ a^^ ^ a^,a.| , . . .8^^. 

Let Oq, ^ the row'veotors of the array g* (S^g}- 
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Xhe noxmal equation (2.2.1) thoi takes the fom 


ir 


3 






( 2 . 2 . 2 ) 


where 0^^ is the (H+1) hgr (M+1) matrix given tgr 



( 1 , 1 ) 




It 18 eas7 't^o see that * (lyl) 
which IB to say 





- 


3 -r,l-l 


•* 


( 1,1 


IPurther, we see that 

If wo define - 

assume^ fom 


A _ A 

r j rt-k, 3 +k* 

0,^ 1 , then the normal equations (2.2.2) 



(2.2.3) 


Xhe matrix appearing In (2.2.3) Is called "block Xoeplits” 
matrix, and we may easily verify that 0^ is a Hermitian matrix. 
In fact, positive definite and k)£o Is Xosplits, but not 

Hexmitiaa in general. 

Our problem is to develop an efficient algontbm fox 
solving the system of equations (2.2.3). Justice [U) has 
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analysed this problem in great detail and he has developed an 
algorithm based on the theory of Ssego polynoDuals. He expands 
the unknown filter sequenoe, a, in terms of the appropriate 
Stego polynomials, to obtain a recursive algorithm for oalou- 
lating a. 


2.3 XtfO SIHMSIQlirAL IXEmrSQir* S ALQOKE!EHH [18] : 

s 

A. Oaloulation of the Bivariate Sjfego Polynomials: 

She ELock-Ioeplits matrix In (2.2.3) consists of 

blocks 

*■ ' Oi3<H, Oik<jr 


She first row of Stego polynomials p^^t ... p^ is obtained 
by the following recursion relation: 


^o,o:o,o " ^^>o 

^o,l+1;o,3 " ^‘*■^^1^ ®o,l;o,l^"^’^^o,l;o,3-1 "^l®o,l;o,l 
where 

X 1 

1 " tio ®o,l;o,k 

x 

Assuming now that the kth row of S4fego polynomials has been 
oaloulated, we oaloulate the (k+1)th row as follows: 


a) 


Cnloulate a polynomial 
^+1,l;i,3 ■ ^,l;i-1fd 


^+1 ,1* 
- Y 

taO 


^kt ^+1 


it;i,d 


- ? 
t^o 


^*kt ^,t;k-l,m-3 


Pi 1< k+1 , 
0 13< M. 
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where 

1 1c+1 N 

^kt “ ♦k+ 1 -r, 3 -B ^+ 1 ,t;r,B^ 

t f k M 

Note -Qiat and depend on the valuea of k and 1* 

Any negative valuea in a auhscnpt of a polynomial ooefflolent 
reault In a zrro value for the ooefficient. 


h) Finally, wo obtain Pj^^^ aa followa: 
^+ 1 ,l;i ,3 " ^k+ 1 ,l;i, 3 /^+ 1 ,! 


where 


f * ^ A 

^+1,1 " '^k+1 tl;k+1 ,1^^Q *k+1-r,l-B ^k+1,l;r,a^' 

Ihia Btep amply noxmaliaea the oorreaponding Ssego polynomial* 

B. levinaon^B algorithm : 

Having generated the orthogonal polynomlalB 

P|^» 04s.k^H, O^l^^H, ae deaonhed in Section 2.31» we 

oaloulate firat iterate o^ to the deaired filter aa followar 

o 

N 

®o.g ■ k ®o,k;o,3 
yhere *ic “ ^ Sq^j ®o,k!o,J 


and la the oroaa-oorr elation array of the input array 

and the desired output array* 
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"b) Assuming kth iterate 
by the equation 

M 

" ®kS>-«3 * Jo 


has been oaloulated to oaloulate 


^kl ^+1,1; If 3 


l^K+1 f 

05. 35.M 


where 


k+l M 

I X 8jpa ^+1,3;r,o 

raO SaO 


®o,3" ®o;o,3 

Ihe desired filter is obtained when Oj^. is obtained* As only 
the polynomials sure needed to oaloulate 0 |^^^ from Oj^f 

the reoursive oaloulation of the polynomials may be done 
simultaneously with the filter ooloulationf the previously 
oaloulated polynomials being discarded as Ihe oaloulation 
proceeds* 

2.4 aianSBAL MOmSL: 

A model which is veiy Important in spectral estimation 
is discussed here. This model is extensively used in spectral 
analysis by Haxiimim entropy method (MEM) discussed in Chapter 4* 
Xhis IS a causal model* Here the two-dimensional signal S 
IS thouc^t to be composed of the message s^j and noise n^^j* 

Xhus we have 

«13 " "13 + “13 <2.4.1) 

where n^^^ describes a discrete stochastic noise process* 
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Qq the other hand the signal is well defined, repeatable 
ginantity* Henoc we oaa extrapolate from the message's past 
history in order to predict with reasonable aocuraoy Its 
future behaviour. 


Thus we can estimate the message at (l,j) from a 
two dimensional fneighbourhood of signals in any quadrant of 
the z-y i0.ane. Xhus 


®J-3 ■ X Jo ® 


'00 


« 0 


( 2 . 4 . 2 ) 


Thus the noise filter for n^^ ^ becomes 
A Q R « 

Jo Jo 


( 2 . 4 . 3 ) 


where 

®m,n " 



^,n 



1 k-l-o 



otherwise 


Let P ■ S[ ^,n ^ noise power. Our aim is now 

to ohoose the ooeffioients He suoh that noise power is 
minimised. She above formulatian again leads to a set of 
equations which eure block Xoeplita. 99ie advantage of this 
formula tion^ called auto-regresslve model (AR) lies In the 
faot that the modol ooefflolente can be suitably dei^oyed In 
determining the power speotral density. 
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Xhe foranilation of the equations and its extension 
to spectral analysis has been discussed in Chapter 4« She 
nuaerioal results of the methods discussed in this chapter 
are given in Chapter 5 » 



OHAPIER ? 


ooirTHncicMrAii angPHQPs op power specikdm ESiimgjiQiir 

Conventional methods for the estimation of power 
spectrum have been in use for quite a long time. In order 
to appreciate the advantages and the short-oominga of some 
of the reoent methods like the Maximum Bitropor Method 
one should have a dear understanding of the conventional 
methods. Ihis Chapter deals with tho two popular techniques 
namely Blackman-Iukey (B-Z) method and the”periodogram” 
approach for the power spoctral estimation. 

She f lj?st step that we enoount' r in th e B-!C method 
IS the estimation of autocorrelation fnnotlon. Ihus in B-I 
method y in order to have a good estimate of the power spectrum 
we must know the autocorrelation values as acourately as 
possible. 

She first section of this chapter deals with the 
estimation of autocorrelation values for two-dimensional 
arrays. She next section discusses the two cooijAentional 
methods for power spectral estimation and their limitations* 
She last section deals with two-dimensional JfBT and also 
two-dimensional windowing techniques. 

3.1 BSIIHallOir OF An!l!OOOREa£U.ZIQir FOnCTIOir (A0F): 

Power spectral densily and AOF are related ty the 
Fourier transform. In some cases even though the estimation 
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of P0«CP .pcot«l densxv xa the B0.1. tlxe 

Of AOP 18 desired as an Inteimediate step. 

Let US consider the estimation of aOr for tso-dimeneional 
data, let I(i,3) the two-dimensional array, let 
denote the tranefl-ated arrays, that is 

I« _(i,3) - jC(i-lc, 3-1) (aero for negative enhsorlpts) 

(3.1.1) 

Ihen the autocorrelation function of the two-dimensional array 
IS defined as 

H^(k,l) - ® a(l,3) Ki-fct 3-1)3 (3.1.2) 

Our prohlen is then to estimate ^7 from a finite observed 
data of siae H by IT. One such estimator is 

0jnr(3c»l) ■ THT i ^ 3[(1-»3) ^ i(i»3) 

** 1-1 3-1 

H H 

■ 1W~ X (3.1.5) 

for aU )k|< M and |l| < V 

She imnediate question that arises lA the goodness of the 
estimator 0^(k,l). A good estimator Should have an axpeoted 
value that is olose to the quantity being estimated. SEhe 
difference between the true value of the ^lantiiy and the 
expected value of its estimator is oalled the bias of the 
estimator. She estimator is said to be unbiased if this 


bias is sero. Trom (3.1.3) the expected value of the 
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estimator Is 

H N 

I [0^(k,i) ]. I i * iKi.j) I(i-k. 3-1) I 

(3.1.4) 

- ^(k.l) for I k l<II , 1 1 l<M 

Ihus we see that 0 (k,l) is an unbiased estimator* imother 

required properly of a good estimator Is that it should have 
a small vici>.inoe. An estimator is said to be oonsi stent 
if it converges in probabilily to the true value of the 
quantity being estimated as the sise of the data becomes 
Increasingly larger, a sufficient condition for an estimator 
to be consistent is that its bias and variance boifti converge 
to sero as M and IT become infinite. 

3.2 POtfKR SEROIRAl ]SS!PIHii!DICII!r; 

Because of the limited observation Interval, auto- 
oorr Nation estimates increasingly deviate from the true 
values as k and 1 increase. Hence k and 1 are normally 
allowed to vaiy between the limits (-K, K) and (-L,L), 
where £ is a small fraction of H and L is a small fraction 
of H. Xhc BLaokman-3!uksy (B— I) method uses the previousHjr 
mentioned estimate of AOB to obtain the power spectral 
density as 

S(u,v) - f f u (k,l) (3.2,1) 

k— K 1— Oi ** 

where j ■ /-I and u and v are the transform variables (e.g. 
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spatinl frequencies in two-dimensional spatial data)* 

A modified speotral estimator 

S'(u,t) k I I 0 (ic.l) Mk.l) (5.2.2) 

kS-I ai-i “ 

IS more acceptable Where hCk,!) is a two-dimensional 'lag 
window' whioh should be properly chosen* S*’(u,v) implicitly 
assumes that AOF is sero outn.de the lag (-fyZ) and (-L^L)* 
Xhis zero extension is artlfl(n.al and degrades the resolution 
performance of the spectral estimator* She two-dimensional 
windows h(lc,l) are discussed in Seotion 3i3b* 


Period ogcom leonnlguft * 

^ Tlu.* vethod emiaqr* the following eatmatoc 

I Jl e-3(»k«l) 1 2 


Ihis method also turns out to be unaooeptable, but various 
windowing techniques improve the estimator* Xhe data window 
h(lc,l) IS chosen to improve stability and mininlse leakage* 
Xhis yields an improved spectral estimator JS^^Uf^)* 


%(u,t) 




I I(k,l) h(Ic,l) 
1^1 


2 

< 3 * 24 ) 


We observe that the term inside the modulus sign is the 
two-dimensional Fourier transform o£ the windowed 

data^ i*e* 
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lOc.l) hdt.l) (3.2.5) 

Xhus 

-inr| * <5-2-6) 

digital oomputer implenentatlon of the Courier transform 
operttlon is best done using Fast Fourier Transform (FFI) 
tcohnlques [7] • Beosuse of inherent periodlo nature of 
FFI, 8 ''jq^(u,v) implicitly assumes a periodic extension 
of the available data. This artificial extension again 
limits the resolution performance. The penodogram 
technique is favoured because it takes lesser oomputational 
time. 


3.3 FAST FOORm TRoITSFQRIf (FFE) jU!n) VUnXWIlTQ IS TSO 
DECTTSIOTO: 

a. Two Dimensional CTI ; 

FFI Is an efficient algorithm for oaloulatlng the 
discrete Fourier transform (BFI) of given data. Ihe UBT of 
a time senes X(n) , n ■ is given by 


B-1 -l*4n 

i(k) * I l(n) 6 " 
n«o 


1c 


(3.3.1) 


Ihe inverse BFI Is given by 

BT 1 3*^^v^** Icn 

x(ii) * 1 1 I'(k) (3.3.2) 

" kM> 
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A dirpot ^pproaoh for oomputatlon of 12ie HBT using (3.3.1) 

2 

requires 4n real multiplications and l!r(4]!r-2) real additions* 

Sinop lAie amount of oomputatlont and thus the oomputatlon 

2 

time IS proportional to IT , It is evident that the nnmher of 
arithmetic operations required to compute the 3)!F9! hy the 
direct method hcoomes veiy large for large values of XT. 

In 1963» Coolej and lEukay 1431 published an algorithm 
for the oompitation of DPI that is applicable When S' is a 
oomposite number* Ihe publication of this paper created a 
flurx7 of aotivily in the application of IVI to signal 
processing and resulted in the disoovciy of a number of 
oomputational algorithms Whidh have oomc to be known as 
Past POurier transform or simpHy PPI algoriiAims. 

Ihc fundamental pnnoip3/P that all these algorithms 
arc based upon is that of decomposing the oamputation of the 
DPI of a sequence of length BB into suooessive small discrete 
Pourier transforms. Xhe manuer in whicAi this principle is 
implemented leads to a variely of different algorithms all 
with oompnmblr improvement in computational speed. Xhe 
prlnoii^e of PPI oan be easily extended to two-dimensions. 

She two-dimensional IPl of an array A(mpn) is defined as 

lUk.i) -“■! V 

miHO n«o 

M— 1 B-xl Trt Irm. 

“"0 1 - 0 , 1 ,..., P -1 


(3.3.3) 



where 


f 




(3»3*3) in a cQ.igh’t*^ difjferenl; wani^aT* 

M-1 N-l n ' w 

I>(lc,l).I IKm.n)!** 1$“ (3.3.4) 

niiBio Ubo " ”■ 


We see that the comput*%tion of 2!(k,l) Involves M jDFI's of 
the form 


A(m,l) 



HmO 


I(m,n) nj“ 


m » 0,1, ...,11-1 

1 * 


(3.3.5) 


followed hy N n7I*s of the form 
H-1 

X»(lc,l)- I A(tB,l) 

TBibO 

1 i* 0,1, ...,8-1 

Ihus the method for eviluating the two dimensional IXfT 
using XPI olgontlm Is as follows. w 


k - 0,1,...,M-1 (3.3.6) 


1. Oaloulate the IVT of individual rows (or oolumns) of 
the given array X(k,l) of else (H ty W) using W point 

\ (or M point) XFI. 

2. Calculate the lOPI of the individual oolumns (or rows) 
of the resultant array A(m,l) (or A(n,k)) using M point 
(or N point) PPI. 

The suhroutine 7ASIPI computes the two dimensional of an 
arr^y X of sise W hy W h7 like above method. lEhis routine makes 
use of nnother routine WS which evaluates the V point SFI, 
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b. Iwo-JDimenBioaal Vindowing t 

Window functions are incorporated Into tbe spectral 
estimator to improve the stability of the estimator. 

Szo^lent baolcground material com be found in BOLaokman and 
Sttkesr [€] , Oponheim and Sohaf or (25] and Eoopmans (W • 

ICbo window w(m,n) should be eiymmetrlo and also have unit 
ampltude at n « 0. Many useful windows have been designed. 
The one used in the pr* sent study Is Idie Seunming window, ietotfa-er 
The ApMndiz gives the routine EAMHIW for this window. Another 
Important reason for using windows is to reduce leakage and 
improve resoluti'n. 8udi a window should have minimum sldelobe 
and niiniim gain lobe width. Both these minima oannqt be 
simultaneously adiieved and henoe a oompromise muet be made. 
Jenkins and Watts [16] point out that the important question 
In empizxoal spootril onolysis Is the ohoioe of window 
bandwidth and not the ohoioe of window function • The ohoioo ^ 

windows in two-dimensions is mu(Ai more oomplioated. Huang (131 
has analysed Ihe proUem of windowing in two dimensions and 
has Uhown that by a suitable transformation, it is possible 
to gft a good two diaenoional window using^good K one- 
dlmenslonal window. 


Let f (x,y) Ibe a two-dimensional array. We would like 
to truncate this array using a window w(z,y). Let 

g(x,y) - f(x,y) w(x,y) (5.5.1b) 

whrre w(x,y) - 0 I x^ + y^l • 
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She Pourier transfoGcm of g(z,7) is 


G(u,v) ■ P(u,v) * W(u,v) (3.3«2b) 

Whore P(u,v) is the Pourier transform of fCz^y) and W(u,t) is 
the Pourier transform of w(z, 7 ). Given a good eTmmstrical 
one dimensional window w*(z)» it Is possible to find a two> 
dimensional window given b 7 

w(z, 7 ) - w»( z^+ 7 ^) (3.3.5b) 


Por ozamplo, if 

v<(*) - M 1*1 

• 0 lz| f. 1 

^u) Is given ’b7 

2 sin u 


W*(u) - 


u 


(3.3.4b) 


(3.3.5b) 


Ihls bfto a first side lobe whose peedc is 0.2 times “ttio pealc at 
u H 0. Ihe oorresponding two-dimensional window using (3.3.5b) 
is 


w(z, 7 ) 


1 

0 


Ix^ + 7^ 

Ix^ + r*i >1 


(3.3.6b) 


!Dhc two-dimensional Pounrr transform has Uie form 


2 2. 


¥(u,v) « 


2ff J^( u +v ) 


u +v 


(3.3.7b> 


whose first side lobe is 0.12 times the petde at u - v - 0. 

A more general comparison would require a oomparison of area 
under the side l&bes in the oase of 1 -dimensional windows'^nd 
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and volume under laie side lolaes in the case of two-dimensional 
windows. Ihe following theorem gives suoli a oomparison. 


Qlheorcm t If a f^ymmetrloal one-dimenslonal window w* Cx) 
and a two-dimensional window w(x, 7 ) arcfrelated hgr 

w(x,y) - w* I (x^ + y^)^l 


then the lE’ourier transforms of w>(x) and w(x,y) satisfy the 
relation 



where 


W(u,v) * H(u,v) 


¥»(u) * H*(u) 


H»(u) - 

H(u,v)- 


1 

0 

1 

0 


U2. 0 

u< 0 

U2. 0 for all v 

u< 0 for all V 


Iho proof of this theorem oon he found in Buang [13] • 


Iho Windowing prooedures are mostly aimed at 
smoothing the spcotral estimate and this causes loos in 
resolution. Ihis is a major problem In the analysis of Aort 
data records. Both B-I and the penodogram extend the data 
In an artlflolol manner without mahing use of any other 
aviillable Information the data provides. Aurg [8] has 
suggested a method of speotrum estimation whioh mahes use 
of the maximuffl entropy principle proposed by Jaynes [15]. 
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IhlB meliiod known as Mamium SitroiBOr Method (MB() extends 
the data In a maxi.tnftlly nonoonmittal mannev outside the 
given data interval. HBl results in a hetter resoluti(» 
and is Idle suhseot of discussion In the next ohapter. 



OHAPm 4 


mnmwi bitobop t mbihdd tMHi) DT jcwD-amoBrsiaas 

As pointed out In the previous ohapter, the oonven- 
tional methods of power speotrom analysis make artificial 
assumptions ahout the data outside the given Interval. 
Moreover, the window chosen for the purpose of Improving 
the stntlstlool properties of the estimate is prefixed and 
does not utilise the information available In the speolfio 
data being analysed. lo overcome these drawbacks. Burg [8] 
proposed a new approach to spectral analysis which has come 
to be known as Maximam Iftitropy Method (MEM). She Barg's 
approach is based on the maxlonua entropy^ principle due to 
Jaynes [15]. 

She Burg* s algorithm has been very profitably 
api^iedbyUlxyoh [31] to solve many problems In geomagaetlan* 
She theoretical considerations concerning the development of 
MM have been presented by Barnard [5] , Edward and 
Eltelson [12], and Emylie et al [29] laooss [20] has 
oompared MM with Maximum likelihood Method. In an Important 
oommunloatlon. Ton der Bos [36] has eihown that ttie extrapola- 
tion of the ACT of a process by maximising the entropy of 
the oorresponding probability distribution fhnotion, is in 
fact eqplvalent to least mean sqpare fitting of auto- 
regressive (AR) model to the process. Using the duality 
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Akalke [3] has given a odtenon for determining the lengtb 
of the prediction error filter. 

In this chapter, we dlsouBs the foxmalation of NBl 
In two dimensions and show how the Wiener filter theoiy In 
two dimensions can he utilised for determining the prediction 
ooefficients. She relationship between the prediction 
coefficients *uid tho power spectrum of the signal is derived. 
S!hc formulation of the^normal equations" and also an algorithm 
for solving these equations has been given. 

4.1 a!WO BIHBrSIOliULL POmSR SFBOXRDH: 

For two dimensional signals, instead of the txme 
variable, 't* ^we usualljr have two space variables x and y. 

In the frequenoj domain, instead of the sincfLe variable *f<, 
we have two spatial frequenogr variables and u and v* Xhe 
oorrosponding power spectrum is denoted by S(u,v). We are 
interested in band limited processes only and hence we define 
A to be the set of all ordered pairs of spatial frequenoies 
(u,v) sudh that {u|£ and |v| ^ V Ihen we impose 
S(u,v) ■ 0 (u,v) }< A (4.1.1) 

!Che power spectrum S(u,v) and Ihe anto-oorrelation function 
P Cx,y) of the oomxAez valued quantity €(x,y) are related 
throufi^ 

+00 

S(u,v) « !! p (x,y) exp {2*i(xa + yv) l.dxdy 

«tOD 

- 1-00 

p (x,y) - // fl(u,v) exp {-2 wi(xa + yv) }aadv 

-00 
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S(UtV) h'^fl the two dimenBional Pouner senes represeat'itlon 
given h7 

+00 

S(u,v)-AxAy I erp{ 2» l(mA X u +^iAy ▼)} 

“.n (4.1.5) 

vhere * P(m Ax, n Ay) are the sampled values of the 

autocorrelation function and 

Ax - 1/2Ujj ; Ay- 1/27jj (4.1.4) 

■ > <«k«5X4« «S,l’ - "X-* 


From (4.1.3) we see that a oomplete knowledge of the power 
spectzum nerds a oomplete knowledge of th e autocorrelation 
values. However, for finite lengtii data, we onn overcome 
this difficultjr by making use of the Wiener filter theoiy and 
the maximum entropy pnnoipae. 


4.2 FOtflSR SPBDIBDN FBOH ADXO-iagGBESSITE HODHL (WISm FIlffiER 
ZHBORI): • 

Let ^ ^ denote the received signal. Ihis signal 

consists of a deterministic message S|g^^ and a random noise 

n||, 2 ^. Xho noise is assumed to be sero mean white. Thus 


•k.l “ "ic,! * "k,l 


(4.2.1) 


Since the message is deterministic an dwell defined, it 
normally {possesses a Tajlox senes expansion in two dimensions. 
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Henop It Is possible to estimate the signal at a point <m,n) 
from its neighbourhood of values In the (x-y) plane, anus 
A Q R 

“ llo llo » *00 " ° 


7or this signal filter, the noise filter becomes 

A ^ A Q R 

\,n " “ ^,n " *hc,l ^m-k,n-l 

Jq ^k.l 5m-k,n-l 


vhcre for n ■ 0,...,Q; 1»0,....R 

k-l«0 
otherwise 




Let P denote the average noise power. 

Then the noise power by emiaoying (4.2.2) is 

2 


^Q+1 ,IUil ■ ® 


? V * I 

Aoio"^ plm ^ in-v1CgZi*J. j 


(4.2.3) 


(4.2.4) 


(4.2.5) 


- I f ^ 2 , 15 w ,, , 1 , ( 4 . 2 . 6 ) 

k,k»-o l,l»-o * K K ,1 1 


Olearly the value of Pq^-| dope nds upon the choice of 
Pk If these are not properly chosen oontalns 

a small faction of the message power also. Our aim is to 
choose I'l. I'B 0uob that the contribution of message power 

Kpl 
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to the aoiee power Is minimum. Por -Hiia «e take the denvmtlve 
®Q+1,H+1 reepeot to and set it sero for nbiO,...Q, 
a ■ 0,«.,pR and m,n both not eero. Ihen wo get 

Itao 1*0 **h-m»l-n (4.2.7) 

Introducing this into (4.2.6), we get 

Q R 

^Q+1,R+1 ■ Pk,l (4.2.8) 

We now define oooffiolents ^ as 

■ ^k,l /^^Q+1,R«.l^* (4.2.9) 

E(],uationB (4.2.6) ond (4.2.7) oan now be written ae 

ItaO llo ‘’k-n.l-o “ *1110 *no (4.2.10) 

mMO,...,Q; n«o,...,R. 


In matrix form this beoomee 
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An algorithm for deteiminlng ijq’b is given In Section 4.5« 
Y* B so dctennined o in be used direotly to determine the 
power spootrum S(u,v)aB shown below. 


4«3 POWXR SEBO'fRUN 8(u,v) IIT lEEHS QFv**S: 

®he noise speotnun ]!r(UfV) Is defined as 


+00 


N(u,v) . 


^ -5) “k«i.l« °i,j} 2 vl(m Ax u-*a Ayr) 

m,n ’ * 

(u,v) e A 


(u,v) 4 A 

(4.31) 


Por uniform white noise iAiis becomes 

, 0 (u,v) i A 

Inserting (4.2.3) In (4.3.1) we have 


(4.32) 


N(u,v) « 


? ! 

k-o L .0 


exp {2 t l(k Axu + 1 Ayn)} 


R 


i I Tici II ®*P{ -2Tr l(k* A « + 1 ' 

IC^hO ^ 


Ayv)} 


Ax Ay I <*Bi,n ®*P 

-GO 

m,n 


(u,v) e A (4.3.3) 


Prom (4.1.3) wo sec that the last factor is nothing but 
S(u,v). IChus using (4.2.9)» (4.3.2) and (4.3.3) we get 
the fln^ result 


S(u,v) 


Ax 4y 


f f ■^,1 Axa+1 dyv)} I 

k»o 1-0 (u^v)e A 


2 

(4.3.4) 
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Ihue "tij obtaining the values of from (4.2.11) and by 
using (4.3.4) the value of the signal power speotxum S(u,v) 
onn be aoourntely determined. 

4.4 SHE VJLUXm BrXROPT 0QHO25Pa!: 

P3?om information theory, the entroxor E, of a disorete 
set, is given by 
M 

S-- I. P 4 I 0 P 1 (4.4.1) 

i»1 * 


where Pj|^ is the prob'lbilily of oooarrenoe of an event m^ 
out of K posoibil tics. Por a oontinuous prooess with 
probnbilily distribution p(q)r this beoomes 


+(r 

H - - ^ p(q) In p(q) dq 

(4.4.2) 

For a Cknassiin prooess 


1 

p(«i) - 7; 7:^ t — 77?* * 

(* v® a r 2 0 

(4.4.3) 

2 2 

whore E( q ] ■ v ; E[ (q- v) 1 •• 0 

(4.4.4) 

She %tropy H then beoomes 


2 ^ 

H - ln(2v 8 a‘) 

(4.4.5) 

entropy fits into the determination of power speotral density 
She problems that arise in oonventional methods of power 
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spootral eBtlmntion ape due to assumptions that are made 
oonoeming the data that lie outside the ohserratlon Intefeal* 
Ihcse assumptions that oonstrain the data to he either penodio 
or soro outside the known interval are incorrect are 
reflected in performance degradation of the spectral estimator* 
What Is required is then an approach that is first of all, 
consistent with the prior knowledge and second, estimates the 
prior probabllltr assignment that describes the prior informa- 
tion without assuming anything beyond that. She relationship 
between entropjr and uncertainty led Jaynes [15] to the 
formulation of maximum entropy principle Which can be stated 
In the following manner: 

"The prior probability assignment that describes Ihe 
available Infoxmatloa but is maximally non-committal with 
regard to the unavailable information is the one with the 
maximum entropy”. 

In order to apply the concept of maximum entropy to 
spectral analysis, we begin wilh the relationship between 
the ontropy And the spectral density of a Gaussian process. 

In all but most pathological oases signals possess a Pourier 
representation. Ve associate to each Fourier oscillator 
expC-2vvt), at frequency v, a oomlex valued Fourier oscillator 
strength f(v). For a stochastic process f(v) is a nndom 
quantify having some probability distribution function. It 
is convenient to employ a Qaussicui distribution in amplitude 
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and n unlfom distribution in phase to independently desoribe 
the osolllator strengths at each frequency* Itaus 

B If(v) ]« 0 ; S(v) « E [f(v) f*(v) ] (4.4.6) 

where S(v) is the power speotrum. Prom (4.4*5)» we have 

H(v) - In I2v e S(v) (4.4.7) 


Wo define the average entroiy over the bandwidth to be 



Given nu to correlation function values ^^+1****'*^H 
the maximum entropy pnnoiple mey be stated as follows: 
Miaximiss Hg with known value of autocorrelation fonotion 
values 


p. ■ J 8(v) oxpC -2x1 ndt v} dv ; n«-M,...,M (4.4»9) 

“ -V 

rur 


She solution to this problem is called the maxliBum entropy 
power spectrum. Since (*^»***> are kaiown we vary 
lnl> N so that Eg is maximisedt that is 

-iSl . 0 p |nl> M (4.4*10) 

»Pn 

Prom (4.4«9) we have 

j_ >LJijydT - 

•Vg S(v) ^n 


0 for ln| > M (4.4*11) 
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Por a band limited apeotium witwaan (-V* and T_) we have 

+00 » a 

I Pn ®*P ( 2 VI n At v) I V |i V- 
nahoo “ • ' H 


S(v) . 


(4.4.12) 


1^1 > ^TS 


IhuB Eg Is an externum when 

+v; 


F 1 

L TT-r e3cp ( 2ir 1 n At v) dv - 0 |ni>H (4.4.13) 


IhuB the Fourier BcrieB for 1/S(v) is txunoated and has the 
foxii 


S(v) 


! 


Bjjj exp(-2v 1 n At v) l"^ 


(4.4.14) 


Sinoo 8(v) In a real quantli^, ooeffiolents are Hexmitian, 
that is, 

■ “in (4.4.15) 


“m 


Epandlng (4.3.4) and oomponng with (4.4.14) we get 

mln(n-4n,m) « 

y ^ 1 jm . . . (M 


(4.4.16) 


Maximum Jntrotar In Two^DlmenBlone : 

We onn extend the maximum entroBjr pxxnolple to 
two dimensions. In two dimensions Eg Is given by 


a iJi— ■ ^ F In lS(u,v) )du dv + ln(2v e)^ (4.4.17) 

®^E^E -t4 -V. 


where 


'N "'H 
+00 


8(u,») . AxAt I P- erp{2»l(.4xa + nAy t)K4.4.16) 

-00 ' 
m^n 
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fti « ■ J / S(u,v) exp {-2ir x(m Acu + iiA7T)>da4!r 

<f ® (4.4.19) 

Haxunm wtronr pnnoiiae than ears that the nhknowi a 

^ BL9& 

Should be varied such that is majdinised. Xhat is, 
d 

-0 |m|>Q or |a|>R (4.4.20) 

Differentiating (4.4.17) we see that Iftie average entropy Is 
in oxtremum when 



exp {2v l(m Ax u -i* n Ay v)} 


S(u,v) 


du dv ■ 0 


|n|> Q or |n |> R 


(4 4.21) 


Thus the Tourlcr transf xm of S(u,v) is truncated and has the 
foru 


S(u,v) - ^ \ 8L. _ exp{ -2v l(m Ax u + n Ay v)> P 

Q n— R » (4.4722) 


Por S(u,v) to be rcil the coefficients a^^jj are Hermitlan. 
Thnt ISy 


«a,n " 



Q ■ '■Qy...fQt n ■ *Rf • . . >R 


(4.4.23) 


The coefficients and are related by the equation 

mln(Q-n, Q) ittin(R-n, R) 


(4.4.24) 
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4.5 AlOORiaaM PQR OAIDUMING WO DDCTSTSIOlttl BIITl® 
OOEPPIOXENTIS: 


IBhc following algorithm oan he employea to solve 
the equation (4.2.11). aChis algorithm is due to Vewmann[24] 
whlhh IB an adoption of Blssanen*s method for solving hlook 
equidiagonal matrioes. Let r^ be the submatrlx given h7 


no 

••• PnR 

« 

• 

• 

• 

• 


ni«>R 

... P no 


n ■ “Q f y +Q 


(4.5.1) 


E- in Hrimltlmt Ihat 1. 
n 



[r*l 


(4.5.2) 


Lrt and be two oolumn veotors of length B+l. 


flk ■ ^^ko' •••» ^IcR^ * 

Ajj, •••! 0]^ 


Ic ■ O^eee^Q 


(4.5.3) 


lEhc equotlon (4.2.11) thou becomes 
-m ■ 1 .^ ^ k-« ^ * m-0,...,Q (4.5.4) 

KwO 


Lot 0® ind rf be two matrioes of size (a+1) by (R+1) given 

O V N 

0° - (4.5.5) 
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ISTow Itomte from N ■ 0, Q-2 

° n -0 ^ 0 nio 

W+1 * -IN 


*^n+1 ^ “ ®n+2J 


(4.5.6) 


0 ‘ gig *n-N-1 ''n ^ 0 ' I -n-1 n -N-2 1 


With for n ■ 0,...,N 


.N+l 

^n+1 




0»+1 

-41 0 


B^'+l - D» - 02 , d"+'' 
n+1 n N"1 0 


(4.5.7) 


0 'ind ]} 'oatPloeB 8‘<tlsi;y 


r 

AbQ 


^K 


®n-k ■ *n+1-k 


N 


i ®k-n ®n ■ *k-n-1 


nao 


k - 0,1,.. .,N 


We dr fine (Rfl) by (R+l) nAtpix P by 


(4.5.8) 


'■‘'o' 5 ! 

IImO 

and D^"'* by 



where I is (N+l ) by (R*-! ) Identily natplx. Ihen we have 

fop k ■ 0».»« »Q 
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Ihc sooond eo^{\tlon is IdmlsLoal 1:o (4«5«4) If we make 1;he 
Asnool'ition 

ftc ■ ^ Ao (4.5.11) 

This algorithm requires 0(Q%l5) operations; to solve (4.2.11). 
Otherwise} wc need 0(W^ R^) operations. 

Iho routine NBM gives Idie Fortran 17 version of the 
'ibovo ilgorlthm. Ihc next ohapter contains many numerical 
fcximplcB nnd discusses the merits and dements of the various 
methods of prediction and spectral estimation disonssed in 
this and the previous Chapters. 
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Bssmas Aim DiscnssipirB 

3!hlB ohaptBP pTBsoats some of 1:lie oompu'tationB.l 
rcoullis obt lined for the xu^ediotion and power epeotral 
cstuQEition nirthode outlined In the previous chapters* 
Section 1 gives the nuucrical results for ABi-HBS model 
for prediction. Section 2 discusses the AlUKEH model 
fitting. Section 3 gives the numerical results for the 
fts well as method of speotrel estimtion. She 
numerical computations were carried out on IM 7044 and 
plots were drawn on 191 1800. 

3.1 AIUNEK KOm OF PBSDIOailOir (SBO.2.4): 

Ilio data us( d for this prediction is the binary 
d)ta of a "girl's picture". She sise of the data is 16 
by 16. 7or this data the prediction ooeffiolnats are 
dcttrmlned using alUMM algorithm. She picture is lAien 
rcc'instructod using these coefficients and the neiebbouring 
points starting with Imown boundary values. Pig. 2 gives 
the. original picture as well as the predicted picture. Xhe 
two levels used in this picture are 2.0 for ^blank* and 
6.0 for 'star'. A better result may be expected for a 
picture with more number of grey levels, as the mimber of 
shirp edges in that case would bo limited. ®ho model order 
used in the prosent case is R « Q " 2. 


Tl-E ORIGINAL PICTURE 
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***•« 
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• * 
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* 
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5*2 MOnEL PimiTG: 

lo test the effioaoj suid oumerloel aoouraq^ of “ttie 
ABf-KM model the scheme as depicted In Iig.3 was adopted. 
Here a two dimensional TNhite noise Is passed through a 
first order difference equation model, with known ooeffi- 
oients (}(k,l), to obtain the two dimensional data z(k,l). 
Uslns the known ooeffloients C^(k,l) the ^tiue speotrum* 
S(u,v) is dctcsmined. !l!he two-dimensional data z(k,l) is 
them passed through the lIUNBf model to obtain the estimated 
parauotcrsG(k,l). Using these estimatee the fHBf spectrum* 

A 

S(u,v) is obtained. She two speotrums are compared. The 
lable 5*1 gives the actual and the estimated parameters for 
various difference equation models of order 1. Fig.4 gives 
the* true spectrum* S(u,v) and Plg.5 gives the HBl spectrum 
^(u,v) for the first model (lEable 3.1), Sl.Vo.l). 

5.3 POWXR SPIOIHQN KBSimDS: 

Xo tost the acouxaGp of the modhl, a two-dimensional 
sincwavo of amplitude 10.0 and centered at frequenegr 
(21^6, 2*^6) was added to white noise of variance 1.0 and 
moan soro. Xhis data was analysed using the above model. 
P],g.6 givoB the spectrum of Ihis data for mod^ order 1. 
Pxg.7 gives the spectrum of the same data for model order 2. 
Xhus we see that increasing the model order gives better 
results. However, there is an optimum value of the model 
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Oiablo 5.1 


Sl.No. 

Oooffloieat 

CKk.i) 

Actual 


H Sexot 

1 

'»00 

0.93989 

0.92401 

3.68 


®ol 

-0.33596 

-0.32321 

3.78 


®1o 

-0.33596 

-0.38784 

14.15 

wm mm < 

0l1 

-0.14398 

-0.13195 

8.33 

2 

®00 

0.93515 

0.87944 

5.9 


0,1 

-0.37406 

-0.36653 

2.2 


Olo 

-0.37406 

-0.32715 

10.9 


0l1 

-0.16833 

-0.12994 

19.8 

3 


0.92420 

0.86025 

7.1 


Oo1 

-0.30499 

-0.30499 

0.0 


Olo 

-0.30499 

-0.26920 

10.6 

a 


Oil 

-0.30499 

-0.30562 

0.16 

■■ ■■ 

4 

®oo 

0.91954 

0.84909 

7 5 


0,1 

-0.59770 

-0.53025 

11.4 


Olo 

-0.59770 

-0.51517 

13.6 


Oil 

0.36782 

&.31170 

15.1 
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order to Iso used* Pig. 8 gives the same slaewave speotrom 
for B-I method. 

Hgsoltttlon Performanoe s 

lo test the resolution perfoxmanoe of B-I as well 

as MBS method, a oomposlte slnewave was used, nhe amiglitude 
of the Slnewave 10.0, !Che two frequencies were (2v/12, 
ar/12) and (2ir/6, 2v/6). She oomposite sinewave was sampOLed 
to obtain a 24 bj 24 data. Ihis was then mixed with white 
noiao of variance and mean sero. The AIU4QM algorithm was 
then used to obtain llie speotrom of this data. It was 
observed that for lower order modele, HBl failed to resolve 
the two peaks (Plg.9» modSl order 1). Is the model order 
was Increased, the two peaks beoame more dear (Pig.10, 
model order 2). Por model order 5, the two peaks were 
dearly resolved (Plg.11)* However, spectral splitting 
was observed In one of llie peaks. Ihls spectral splitting 
has been observed by others also (44). Pig.12 gives the 
same composite sinewave spectrum for B-® method. 




tier 


Pig. 4s frue Spectrum of the signal x(k,l) generated "by the 
first order difference equation model (using the 
true coefficients G(k,l)), 







lig.S Spectrim of sine wave centered at frequency [2s /6f 2s/63 
using B-T aetBod with Hanning window# 





a^olutloD Of oo»5,oslte sine save of tso fteque 
order* 2*/l2] losing ae-mem model 




Resolution of composite sine wave of two frequencies 
[2n/5, 2n/5] and [2%/l3t 215/133 lasing IR-MM model of 


order 5. 


Pig. 11 




Hg.l2 Resol-ution ©f coaiposite sine ware of freq-uencies 

[2%/6, 27t/6] and [2^12, 2icA2] rising B-f method 
with Hamming window. 
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APPENDIX A 


This contains the Wiener filter synthesis, hvased 

on theWiener filter theory discussed in Sec. 2.2. The 

algorithm used here is the one given in Section 2.3. 

This algorithm could not he implemented for solving the 

prediction equation (4.2.11) as the term comes on 

the right hand side and we cannot evaluate Y directly. 

00 

However, two examples of Wiener filter synthesis, using 
this algorithm are given helow. 


Example 1 

Input array to he filtered: 

I 1 4 -1 

h = j 

I 4 0 6 

i 

Desired output array: 

1 -1 ~3 1 

4 _2 10-8 

c = 

0 15 7 11 

12 4 13 6 

The autocorrelation of input array h: 


58 

0 

-2 

12 

0 

0 

0 

58 

8 

-2 

0 

0 

-2 

3 

58 

0 

-2 

12 

12 

-2 

0 

53 

3 

-2 

’ 0 

0 

-2 

8 

58 

0 

0 

0 

12 

-2 

0 

53 
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The cross-correlation array G is given hy 


i 32 

-62 

32 

152 

179 

170 


The synthesised filter array: 


/V 

a = 


1.015 -1.031 

0.595 2.023 


2. 327 1.153 

The actual filter -- ’'a*: 


a 


ri 

s 

3 


-1 

2 

1 


Example 2; 


Input array: 


h = 



2 3 

5 ‘ 6 


Desired output array: 


j' 1 5 9 9 

3 17 24 24 

-2 312 15 

8 4 17 6 

- 
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Auto-correlation matrix: 


<P = 


I 91 

I 

53 

32 

17 

0 

0 

i 58 

91 

83 

32 

0 

0 

32 

23 

91 

58 

32 

17 

I 17 

J 

\ 

32 

58 

91 

23 

32 

! 0 

0 

32 

23 

21 

58 

I 0 

0 

17 

32 

58 

91 


Gross-correlation array: 


G = 


267 

207 

254 


372 

315 

254 


Sunthesised filter array: 

r T 

: 1.0 3.0 j 

•a = --I .O 2.0 1 

I 1 

; 2.0 1.01 

The actual filter: 

1 3 

a i -1 2 

2 1 i 



APPENDIX B 


Phis contains some of the important listings of 
the computer subroutines used in the computer simulation. 
They are WHITE, C07AR, WUNIE, SPBOTR, EASTPT, PET, H&MMIN 
and MEM. Comment cards introduced in each routine explains 
the function of that routine. 



c 

c 


$0 

lOo 


IHIS SMBRCUTlHe <s€N£RATiS W ‘rjp: cf Vp 

N SV N-MUH A KNMH «6A« ANB VARiAStsi 


N BY N-Ut7fi M: m(mi4 mAM Am- VAR'Ijt’r^'ti* 

0 1 mm 1 ON X f 62 s i » w i 2 s , 25 r ' , 

kIkIT^SoO 

UMVAfi«VARWS/i2, ' ■ 

H=fi2*N)+2 

C^IL wyNlFCy,H^UW|V.AR,RIfs{|T) ■ '■ 

L*l * , ' " / 

OC l-:iM K=sl,N 
X I K ) s y a i,< 

Le=i+ji . 

OC uJ i»L,LS 
X(KNX{K)+YC £) 

- Lal + 1 ;> 




, ^ 00 30 ]=]Is 25 

• -00 30 J«l,25 

K- 1 1 *11 -4=2544!, 

30 WCi,J)«X(K} 

■'.RETURN ■. 

END’ ■ 



SUBROUTINE CCVM(A,C,HArfe>t : ■. 

THIS SUf^RGUTi-Ne CAE^^tgS. T««‘ 

Array X OF SI2e H FOR SWIFTS.' 

THE AuTO-CORi^tltt.l-g^^ VAIUFS'-ArE -Sf6ktj-4’#;; 


''G, 4? 




'i ^41*1- C, t S , 3 1 ■, u C 24» 24'K 


Ci r> o o r> 


CALCUtATiNG THE SPfCTRAl VAtUE^S S FROM 
L CUEFFICIFNTS G.-Un»VN DENOTE THE .NUMBER OF POINTS 

•.M iHr SPECTRAL GRID* R,Q DENOTE THE ORDER qF THE AR-MtM MOoeL* 

DIf.»EME |n.'4 S( 32 , ,G(2i 2) ' : 

1 M EGtR Rf Rl,Ry ^ Q, Ql,g2,u?4tVN 
CGf'PL5=X SUMtZ.XX 

P i ■ 4 j " ’ 

R 1 =R+ I 
Ql.sO+i 
R2 = 3*R+1, 

OElX=i^/{2BAFLCAT{UN)) 

DFLY=|*/( 2 , «iFL jAT{ vN ) ) 

OC 1 ilJU=J,UN 
OC 1 iVV=J.»VK 
lU^-lUU-l 
!V«JVV-l 
SUM«{w).,0*} '-v 

DO 2 KKalf 

nu 2 LL = l,tQl :. 

K®KK-1 

L-LL-I. , ; 

YY»«pLOATCK )*CEl.X*FtCATC iUI + FLOATCL }1'OeLY»FLOATfI vr 

Y«- 2 * Yf-’i^YY • ' 

XXacMPLXtO, ,Y) ■' ' ■ "" 

z»cexp(xx) 

SUP»SUM+Gf KKf LD -^Z 

SUUU, IVV)aC0ELX<'0ELY)/{CABS{SUM)**2J 

RETURN' 

END 



SUeROUT'tNi: NEMU,R',0,G J, ' • I,/ " 

SUBROUTINE FOR CAiCULATING TMI' COtFFlClBJ^j'' ■■f3N^>HE':'A|^|«eM‘ 
A IS THE INPUT AUT O-COftREL ATIOM 'MATRI xi' R^ANo ' q OEI^ffi THI 
OF «yT^^ 


4 . 

'-*=■ '-'C + O 

DC “ -Ks-URl 

OC L^'V,ri 

0 C rs Kf L ) -*3 <, 

C ^ P f K f L ) X 0 « 

5 

CCKTiOUC 

20 

CUTAt 7 )«i 

OC ?0 

OtLT/ (1 jr^: 

30 

OC 2u 

OC 10 

Alf'iVf If J) = ACC1 


Ci-il l••ATi /il . Vji-j , e, 3, 0 = 1 

p s 


DC 4'T 1 = 1 , RI 
Of* ^'0 

K»l, +Q-i-J 

^0 DUiwi<l,J)x^{K, || J) 

DC 4« I = t,Rl 
DC 4'5 J'-ijfti 

^5 . nu^'2(l,, = . 

CALL f^MLTH OyM;!t r;UMl,Nf H) 
OC hO i=|,Rl 
DC S»0 

5'0 C(P,1,4}«0UM2C i, J} 

OQ 60 IwljRI. 

DO 60 J*UR3. 

K»C 



DUWiatO)*AiKf StOJ 
DO 65 . ■ 

DC 65 J«lf R| 

DUP2f |,J)«AINvn,J ) 

CALL MULTI (OUM2rOUH1,NtM) 
DO 7a ,!:,-URJ 
OG TO •' 

0{P,I,J)«0UM2Ci,J) ■' 
iPco.Eaa) GO TO ?? 

OC 80 NNl»i,G3 

NiaNNi-i- ; 

Oi*»NNi+I , - ^ 

00 85 'I-i, RI 
00-85 

dum 3 < f f j '■ • 

00- '90' HW2®it 01": ■ 

P*N 2 *Q*H»i:. ••■■ 

KsKNI-N.'+Ql 
OC 

m: 

DUMU IrOlaAIK,'!. 

ouMoc 

CALL PtllLTl'mii 
00 9*^ i-^um 

0 L 9 0 •0”* X t ■Ri| 

CKif^M i».o 
CawTii 4 MC 


11 0 

115 

130 

125 

120 


ivw 

rjuwj (:;,JU-HKy i,J) 

L«'j 

CALL A0D(C!JM^,DU^*Ut,^U 
C41L f^ATf'^Vf CyM3?M»8f0,DET?RMI 
Ilf l-URi 
DC ■■ ■' •S J=:1,Rt 

on t"ij 

Nisn'-i;-,'-’ 

K-^^zh'^y+Q 

3'.- 13 1 

DC j 3 I J-i,R'; 
0UHltj,}J)-,3{XtifJ) 

0U^»2Uf J)»CC3, V, j) . 

CAIL VUl'fICDUHi,0UH2»M»H) 

00 12 5 i:=.-l,Rl 
00 123 J = 

Dy)«-!t( % j) = dum3CX9J )40um'.( i,j ) 

CCMIMJE 

KaQl.+KNi4i - 

OG^UO 

DC 140 J»Uri 
OUMit |» J) = ACKf I, J) 

CAlL, MATAODC CUM4fDUMl,L»M) 

CALL MULT|«DUMi,0UM4fH,Mi 

PaKNI 41 



DO 150 i»l,Ri 
00 150 J«l»Rl ' 

C{pf?;» J)«f>UM3C It 03 
00 19:5 f = l,Ri 
DC 155 JssifRi 
0UM3( If J }='(!•« 

Oq lAO NN2«l*J?, 

N2«NN2-1 

P«N2*0+flwi 

K*N2-NNJ.4G1 
OQ 170 

00 }.To J=i.,Rl 

OUHl{ K, 1,03 ■■■' ., 

CALL 

ft: '.?«• ; = 3»R? 


1, f* 


Hi 

, j'j i i s j 1 3, : 

■)': lo:. /‘;n 3=; jH-5. , ' , . 

fj.i.j’s 41, 

rr 7j)o 

nC' ji)i} J “ I. f H I 
0U^<^( I J J) = 4( K, i, J) 

20u 'OUK^C , f Ja«C{ f», it J) 

CALL WULTICDUMi,,DUM?fNfM} 

OG 205 I»3.,Rl 
00 Is1 

205 0U^«4( If, J)«DU«4(ItJ HOlWli i, J ) 

190 CONTI MUe ' ' . ' 

K«0l-(NN!. + 13 
00 210 l»l,R,S „ 

DO 21t»-' 

210 DUKlCIf I,J) 

CALL MATA00(BUH4,0UMltLf?<!) , 
CALL f^ULTI(OUM.1tOUH4#.NtH} 


P»NNl+l 



00 220 

oC'P,#-i»j}«oumu» j) 
O'O 230 

P«NH2>»Q-»-NNl'»'i 
Pl-NaTQ-t-NNl 
P2«(Nl“fi2)*0*NNl 
P3*NNi + i , . ■ 

00 24p tHk.Rl 
00 240 

DUHK If ) 

0UM2Ut J)*C(PBfl'tJ } 

I 


CALL' MULT|TOljM5,OUH2»M,M) 
-'CALL IfA^rAO'CCCuMOtDUMlfLfW) 

■W-m4 

.;o-o 

mQ. 240- 






r: 

/. F 1 FQ^.J ) GC /{ 0 225 
'F I, F s i f J 1 =*’< <3 



■o j n 24 3 

2 53 

L lF»T}J}=— j, 

243 

CONTI i\ Ur 


DC 265 l»l. 


00 265 J=lf 

265 

'0UM3C If JJsU 


DO '27 3 KN1«1«0 

p«!^l«C+Q.- 





320 


Kao-*.ji+ai , ■ , "... 

DC 280 ■ 

DC 280'.=J*l,ITi . ■ . . 

ntJHiCI, JJaAlKf h j) 

CALL MULTI C CUMI f aiM2fNfH) ■ 

Ou 29A |=.1,R1 
DO Z.9() J»i, Ri' 

0UH3C I, = ) + OUMl( I, J ) 

CONTINUE , . " ■ 

DO 3tK* 1*1 fRl 
DO 300'J«lyRi ' 

DUMlC If |» J) 

CALL ^ MAT, ACD't€UMlfOUM3fL#N) ' 
CALL Mat INVC OUMi,^M, BfO, DET fiRM I 

DC ?>iD I!*ltCI . 

" 

OC 320, R*lfR| ' 

DQ 320 L«l»R| ^ 

■ 0UM2(K,U*CCPtK,i) 

CALL MULTI<,CUH2,DUM1,H,M) 

CALL VECTOHC CUM2f0gLTA#0UT,NfM) 
00,33U ■ 

G«I2,K}*-#UT{K) ■ ■■ ' 

, CGf^TTWe''' ^ . 

f , ' t* F , '' T 4 






SUBftOyTW^' FFTC-ArNtHU) 


-.UtiPv'Ur !NU calculates 1H:. 2-0 FFl Of ThE' input array A:, 

C.CLPLCX AINWtf|V»),j>(32),Q<'i. JsTEMPOi) 

M = 5 ■ , ■ 

UL .litjy i~.i M . 

DC K=i-tN 

TCpp(K)“AC i,K) 

call FFTnCMPyf^jM) 

OC yjyZ 

AUsK)=TeMPIK) 

CONTINUE'" 

DO 3nu "y»ltN ■' ' 

0C:^30l 

TBMPfK=f«A-|Kf4) ■ 

..■XALL-’FFTCTEHPfN*H) 

■ 0Q.30|.^ ^ Dft' ' ' ^ 

AlK#41*r«WtKJ '-C ■ 
continue > 

RE'TURN ' • ■ :' 

END ■ " 'A h;'- 


THIS' SUBWOMTiNt COMPUTES ‘the F?T XF A " si6UE»lCr"’A OF LENGTH N*‘ 

^ -■■■' X^EMENGTH N SMOULC BE A ■ POWER Qf-i, NO IS LOGARITHM OF N TO; I 


OkFgNSION XREAK .lOG ) jX i?-5AC( 100) 


;0'0 10 I»t,'N ' ■ ' 

■'■:D ■KWGlIJsAIMAGlAa)) 
''■«R£AL(ir*R|ALiA« I) ) 
C0^T!^4Ue D" : , 

• N2-N/?. . r- 

NUl«NU-l 

Kwi’l •' T ' ^ 

■ OQl0Ol=lf NU ' '' 

001:01 l*UNa- 
P»1.8ITRIK/Z*fi 



p ( .If GnTrr-. 

< j 

-1 

lOj 

r,P'i'r^K =:.5 N 

;-i& r.RfK-;, NU) + i 

ifii^LF.,K}' GOTO la i 

T,«E.it=XHFALCK) 

n.MAG*X!,MAGCK} ■ 

Xa64LCK|*XReALf|l' 

^ XIMAGCK)aX|«AGiI J 
, XHEAi(|J*ftEAt ' " 
XI»AG«n«fiMAG.,. '■ 
lot ... COWTINUE:;. ■ ../*• ■ ' 

‘ ■ R'l*XReAl<-|;i 
Ai«xiMA6ii,i'"::v 
104 A(I 

RETURN ■ - 

6 N 0 ' 


V ’ 4 ■ 'V'.’l.,. • 

‘I- 

42 «' 4 I/^ 

■ ' : lBlt,R»lSiTR42-“»-i41-24j'2^ 

200 4 |w 42 '.''''<!'n>flfe.::::''‘ 

. , RETURN ' '• ’■■a; I'/. 

, : ENO ' . ‘ , I. 



